Comparison of low amplitude oscillatory shear in experimental and computational 

studies of model foams 
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A fundamental difference between fluids and solids is their response to applied shear. Solids pos- 
sess static shear moduli, while fluids do not. Complex fluids such as foams display an intermediate 
response to shear with nontrivial frequency-dependent shear moduli. In this manuscript, we conduct 
coordinated experiments and numerical simulations of model foams subjected to boundary-driven 
oscillatory, planar shear. Our studies are performed on bubble rafts (experiments) and the bubble 
model (simulations) in 2D. We focus on the low-amplitude flow regime in which Tl bubble rear- 
rangement events do not occur, yet the system transitions from solid- to liquid-like behavior as the 
driving frequency is increased. In both simulations and experiments, we observe two distinct flow 
regimes. At low frequencies ui, the velocity profile of the bubbles increases linearly with distance 
from the stationary wall, and there is a nonzero total phase shift between the moving boundary 
and interior bubbles. In this frequency regime, the total phase shift scales as a power-law A ~ oj" 
with n « 3. In contrast, for frequencies above a crossover frequency uj > ujp, the total phase shift 
A scales linearly with the driving frequency. At even higher frequencies above a characteristic fre- 
quency Uni > Up, the velocity profile changes from linear to nonlinear. We fully characterize this 
transition from solid- to liquid-like fiow behavior in both the simulations and experiments, and find 
qualitative and quantitative agreement for the characteristic frequencies. 

PACS numbers: 05.20.Gg,05.70.Ln,83.80.Iz 



I. INTRODUCTION 

Aqueous foams are collections of gas bubbles that are 
separated by liquid walls [l|, and like other complex 
fluids, such as pastes, emulsions, and granular media, 
they exhibit transitions from solid- to liquid-like behav- 
ior in the response to applied stress or strain. For small 
strains, foams behave elastically with stress proportional 
to strain. Above the yield strain or stress, bubble re- 
arrangements occur and the system behaves as a liquid. 
In contrast to Newtonian fluids, foams display complex 
spatio-temporal behavior in response to applied shear in- 
cluding intermittency, shear banding, and nonlinear ve- 
locity profiles 0, i, i, i, i, 0, i, 1 E [O, m • Despite a 
number of experimental, numerical, and theoretical stud- 
ies of driven foams, a fundamental understanding of the 
response of foam to applied shear is still lacking. 

In this article, we describe a coordinated set of experi- 
mental and numerical studies of model 2D foams under- 
going applied oscillatory planar shear to characterize the 
transition from solid- to liquid-like behavior and from lin- 
ear to nonlinear velocity profiles. There have been several 
studies of the response of foam to steady shear, however, 
most of these have been performed in the Couette ge- 
ometry in which flow is conflncd between two concentric 
cylinders 'l3|. Instead, we focus on planar shear to avoid 
the 'trivial' transition to nonlinear velocity profiles that 
stems from the fact that in the Couette geometry the 
shear stress varies with distance from the center of the 



shearing cell. 

Another distinguishing feature of this work is its fo- 
cus on oscillatory rather than steady shear as the driv- 
ing mechanism. There are several reasons for selecting 
oscillatory shear. First, oscillatory shear allows one to 
control the amplitude independently from the frequency 
of the driving. When foams (and other complex fluids) 
are driven by steady shear, they exist in a highly flu- 
idized state that is characterized by continuous, often 
highly correlated bubble rearrangements, or Tl events 
[lji|. In the highly fluidized state, the statistics of Tl 
events determine the flow curve and control stress fluc- 
tuations [llIiEHEHIillio, 21]. With oscillatory 
shear, one can study the low-amplitude flow regime in 
which particle rearrangement events do not occur, yet 
the system can transition from solid- to liquid-like be- 
havior and from linear to nonlinear velocity profiles as 
the driving frequency is increased. Since Tl events can be 
suppressed when using oscillatory shear at low amplitude, 
the dissipation between fluid films becomes the dominant 
relaxation mechanism |22l |. Thus, in this regime one can 
directly probe the dissipation mechanism by tuning the 
driving frequency. 

In this article, we report on combined experiments and 
simulations on model 2D foams: bubble rafts in exper- 
iments [23j and the bubble model in simulations [24|. 
Bubble rafts consist of a single layer of bubbles floating 
on the surface of water. (See Fig. [T] for a snapshot of the 
bubble raft used in experiments.) Bubble rafts have a sto- 



ried history as 2D model systems for both crystaUine and 
disordered sohds [23,[23]- In addition, we have performed 
a number of studies characterizing these model systems 
by measuring and quantifying Tl events 261, stress fluc- 
tuations [27|, velocity profiles [a, S S 3., 23., 22] , and flow 
transitions [8|. In this work, experiments on 2D bub- 
ble rafts will be compared to simulations of the 2D bub- 
ble model introduced by Durian [24]. The bubble model 
treats foams as soft disks that experience two pairwise 
forces when they overlap: a repulsive linear spring force 
proportional to bubble overlap and a dissipative force 
proportional to velocity differences between bubbles. A 
useful feature of the bubble model is that it can be gen- 
eralized to particles with finite mass [y, [D]. Thus, the 
ratio of the damping and inertial forces can be varied to 
interrogate the damping mechanism. Recent work has 
shown that the bubble model successfully captures some 
of the features of the dynamics of bubble rafts under 
shear, for example, the statistics of individual Tl events 
[26|. Thus, a comparison of experiments of bubble rafts 
and simulations of the bubble model in a well-controlled 
planar shear geometry will allow us to further test under 
what conditions the bubble model accurately captures 
the dynamics of model foams. 

We will focus on measurements of the total phase shift 
between the driving wall and interior bubble displace- 
ments, and velocity profiles in systems subjected to low- 
amplitude oscillatory planar shear. At low driving fre- 
quencies w, we observe a non-zero total phase shift, while 
the velocity profiles rise linearly with distance from the 
stationary wall. At low frequencies, the total phase shift 
scales as a power-law A ^ lo^ with n « 3. In contrast, 
for frequencies above a crossover frequency w > Wp, the 
total phase shift A scales linearly with the driving fre- 
quency. At even higher driving frequencies coni > ojp, 
the velocity profiles transition from linear to nonlinear. 
We compare the two crossover frequencies LOp and LUni in 
the experiments and simulations and find both qualita- 
tive and quantitative agreement. The structure of the 
remainder of the manuscript will be organized as follows: 
section II, theoretical background; section III, simulation 
methods; section IV, experimental methods; section V, 
experimental and simulation results; and section VI, con- 
clusions. 



II. THEORETICAL BACKGROUND 

We will now review a simple theoretical treatment of 
the response of an idealized viscoelastic material to an ap- 
plied oscillatory strain, which will provide a framework 
in which to interpret the experimental and simulation 
results in Sec. |Vl The main point of this section is to 
identify the possible flow regimes in viscoelastic materi- 
als and their distinguishing properties. For illustration 
purposes, we have selected the Kelvin- Voigt linear vis- 
coelastic model with frequency-independent elastic mod- 
ulus and viscosity, though we will discuss how the results 
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FIG. 1: A typical snapshot of the 2D bubble raft used in 
experiments with packing fraction <j> ~ 0.86 and N m 3700 
bubbles. The bubble raft is composed of a bidisperse dis- 
tribution of bubble sizes: a 4 to 1 ratio of 2.5 ± 0.3 mm to 
5.3 ± 0.5 mm diameter bubbles. 



from this model can be generalized. 

If a planar oscillatory shear strain is applied to a vis- 
coelastic material (with shear along x and shear gradient 
along y) , the x-displacement field Ux (y) of the system rel- 
ative to the initial positions can be obtained by solving 
the force balance equation [30| : 



as 



xy 



dy 






(1) 



where the velocity field is v^ = dux/dt and p is the areal 
mass density. The shear stress 'Sxy includes both the 
elastic and viscous contributions. As mentioned, we will 
focus on a linear viscoelastic material with 



G-f + T]j, 



(2) 



where the elastic contribution is proportional to the shear 
strain 7 = dux/dy and the viscous contribution is pro- 
portional to the shear rate 7 = dvx/dy. G is the elastic 
modulus, and 77 is the dynamic viscosity of the material. 
In general, complex fluids possess complex shear moduli 
with arbitrary frequency dependence. We have consid- 
ered this case, and find that the quantitative scaling of 
the crossover frequency Up depends on the details of the 
viscoelastic model. However, the qualitative features of 
the Kelvin- Voigt model, i.e. the existence of LOp and w„i, 
are robust. 

We consider the case of parallel plates aligned along 
the X-axis separated by a distance Ly in the y-direction 
as depicted in Fig. [TJ The boundary at y — Q is station- 
ary Ma;(0, t) = (bottom boundary), and the boundary at 
y = Ly moves according to Xb{t) = Ux{Ly,t) = ^sin(wt) 
(top boundary). The same geometry and notation is used 
for the experiments and simulations. To solve Eq. [I] we 



use the ansatz Ux{y,t) = A{y)sm{ujt) for the displace- 
ment field. Putting these elements together, we find the 
following solution to Eq. [1] for the displacement field 



and 



Ux{y,t) = Im 



v.j:{y,t) = Im 



Ae 



tu^t sin(fcy) 
sin(fcLy) 



Aiuje^' 



siTi(ky) 



sin{kLy) 



(3) 



(4) 



for the velocity field. In Eqs.[3]and|4l the wavenumber k 
is complex, and satisfies the dispersion relation 



Gfc2 



U! 



,rjujk'^ 



^V~p) 



(G2 + (77a;)2)i/2 



(5) 



(6) 



Distinctive features of the velocity profile are best de- 
scribed by rewriting Eq. |4]in terms of a j/-dependent am- 
plitude Vjnagiy) and local phase S{y): 



Vx(y,t) ^ V.mag{y) COs{ujt - S{y)). 



(7) 



We define the total phase shift A = S{0) — S{Ly). Because 
the flow is periodic, the velocity profile at a given time 
t repeats at subsequent times separated by period T = 
2tt/lu. To simplify the analysis, we will focus below on 
velocity profiles at times when the boundary velocity is 
maximum (i = and Vx(y,0) = Vmag{y)cos{S{y))). In 
the simulations and experiments, statistical accuracy was 
improved by averaging over driving cycles. We defined 
Vx{y) = {vx{y,2Trp/uj))p, where (.)p indicates an average 
over p cycles. Monitoring the full time-dependence of 
the velocity profile is important, but is outside the scope 
of the present work. Error bars on the local phase shift 
and velocity profile in the simulations and experiments 
are given by the rms fluctuations within each bin and are 
typically the size of the data points in the figures unless 
otherwise noted. 

It is instructive to consider two limits of the dispersion 
relation in Eq. [S] the limit of a pure solid {G ^ 0, i] — 0) 
and the linut of a pure liquid (G — 0, rj ^ 0). For 
the pure solid, we recover the dispersion relation uj/k — 
\J G I p = Vs , where Vs is the speed of shear waves in the 
solid. In this case, k — lo/vs is real, and the velocity field 
is a standing wave given by 



Vx{y,t) — Aujcos{ijjt) — 



sm{ujy/vs) 



sm{LuLy/vs 



(8) 



For w < uj'^^i = Vs/Ly, sm{ujy/vs) ~ ^y/vs, and 
the velocity profile becomes linear in y/Ly, Vx{y,t) ~ 
Aiv cos{u!t)y/Ly. For fixed system size, the transition 
from linear to nonlinear velocity profiles occurs when 
u! > oj^;. Because k is real for the case of the pure solid, 




FIG. 2: (Color online) Normalized horizontal velocity profiles 
Vx(y)/vx{Ly) at i = for the (a) pure solid and (b) pure 
liquid obtained from solutions to Eq. [T] as a function of the 
driving frequency u. We show ui/uni = 0.1 (squares), 1.1. 
(circles), 2.1 (upward triangles) , 4.1 (downward triangles) , 6.1 
(diamonds, only in (b)), and 20.1 (left triangles, only in (b)). 
When referring to the solid (liquid), ujni corresponds to uj'^i 



the total phase shift A == 0, and the system oscillates in 
phase with the driving wall. 

In the limit of the pure liquid, we recover the dispersion 
relation iuj = —vk^, where v = ij/p is the kinematic vis- 
cosity. In this case, fc = (1 — i)/D with D = ^ 2r\ j (lo p^ . 
The form of the velocity profile is more complex than 
that for the pure solid. However, for small driving fre- 
quencies the velocity profile can be expanded in powers 
of Ly/D, and the first term is linear in y/Ly. Thus, for 
Ly/D < 1 or w < 2w^;, where w^; = ri/{pLy), the ve- 
locity profiles are approximately linear, as we found for 
the pure solid. However, in contrast to the pure solid, 
there is a non-zero total phase shift A in the pure liquid 
since the wavenumber k is complex. The full form of the 
phase shift is complicated, but at low driving frequencies 
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FIG. 3: (Color online) Normalized horizontal velocity profiles 
Vx{y)/vx{Ly) at f = for a viscoelastic material with /3 = 
GpL^/rf — 1 obtained from solutions to Eq.[T]as a function of 
the driving frequency uj. We show u/LUni = 0.1 (squares), 1.1. 
(circles), 2.1 (upward triangles), 4.1 (downward triangles) , 6.1 
(diamonds), 8 (left triangles) and 20 (right triangles). uOni 
corresponds to o;^; — uj'^i. 



one can expand A about a; = 0. For the pure liquid, at 
lowest order in a;, we find that the total phase shift scales 
linearly with the frequency, A ^ ijj/Quj\^i. 

In Fig.[2l we plot the velocity profiles that satisfy Eq.[T] 
for two cases: (a) the pure solid (G > and -q — Q) and 
(b) the pure liquid (G = and 77 > 0). In Fig. [H^a), we 
show that the velocity profiles for the pure solid become 
nonlinear when oj/iOni > 1- Note that the profiles first be- 
come nonlinear by developing negative curvature above 
the linear profile, and then as the frequency is increased 
further the profile develops positive curvature below the 
linear profile. This nonmonotonic behavior is caused by 
the standing wave solution in Eq. [51 In Fig. [Hb), the 
velocity profile for the pure liquid begins to deviate from 
a linear profile when oj > 2uj1^i. In contrast to the pure 
solid, the liquid system only exhibits monotonically de- 
caying nonlinear velocity profiles with a 'decay length' 
that decreases continuously with increasing driving fre- 
quency. An interesting feature of these profiles is that 
one can observe negative velocities at sufficiently high 
frequencies as in the case of the pure solid. 

We now consider the solution to Eq. [T] for the velocity 
profile in the more general case of a viscoelastic material 
with nonzero G and 77. In Fig. [3l we show an expanded 
range of driving frequencies for a viscoelastic material 
with /3 ^ «,/c.i;)2 = GpLl/Tj^ = 1, so that 0;^; = w^,. 



Similar to the case of the pure solid, when lu 



> 



''nl 



the ve- 



locity profiles show a small negative curvature with the 
profile shghtly above the linear profile at w/uJ^i = 0.1. 
At higher frequencies, the system behaves similar to the 
pure liquid with monotonically decaying profiles and a 



FIG. 4: (Color online) The total phase shift A = 5(0) - S{Ly) 
for viscoelastic materials with different values of /3 = GpLi^/rf 
versus the normalized driving frequency uj/lo\i. The curves 
for all /3 scale linearly with u at sufficiently high frequencies. 
We include /3 = (solid line), 0.0005 (squares), 0.005 (upward 
triangles), 0.05 (downward triangles), and 0.5 (filled circles). 
The crossover frequency uOp can be obtained by locating the 
intersection of the low-frequency and high-frequency power- 
law scaling forms for A. The slope of the black solid (dashed) 
line is 1 (3). 



continuous decrease in the decay length with increasing 
frequency. At sufficiently high frequencies, we also ob- 
serve a regime in which the velocity becomes negative. 
For viscoelastic materials, we define ujd > ^ni as the fre- 
quency above which the system begins to display mono- 
tonically decaying 'liquid-like' velocity profiles. 

In Fig. m we show the total phase shift A for the pure 
liquid and viscoelastic materials as a function of the driv- 
ing frequency. (The pure solid is not shown since A is 
identically zero for all frequencies.) As expected, A for 
the pure liquid (/3 = 0) scales linearly with uj. For vis- 
coelastic materials with /3 > 0, there is a clear crossover 
from low-frequency scaling A ^ uj^ to high-frequency 
scaling A ^ lo"^ with m < n. For the Kelvin- Voigt 
model, the high frequency limit is equivalent to G = 0, 
and A — oj/6llj1^i, corresponding to to = 1, and the 
low frequency limit is A = a;'^a;'j^;/(6wj^;*), correspond- 
ing to n = 3. Using these expressions, one can derive the 
crossover frequency explicitly, ujp = w^j^/w^; = uj^iP- 
For a more general model with a complex shear modu- 
lus, where the stress is given by Hxy = G*{uj)^, G*{uj) — 
G'{lu) + iG"{uj), and G' and G" are the storage and loss 
moduli, the values of n and to depend on the frequency 
dependence of G' and G". However, for physically moti- 
vated G*{uj), the crossover from the low-frequency elas- 
tically dominated to high-frequency viscously dominated 
behavior persists. One consequence of the crossover in 



frequency dependence is that the low-frequency total 
phase shift tends to zero rapidly at low frequencies, and 
thus A may be difficult to measure at low frequencies in 
experiments. In our experimental studies, we were able to 
detect the change in scaling behavior of A, but were not 
able to measure the scaling exponents accurately. Much 
more sensitive experiments are planned to measure LUp 
and the storage and loss moduli at low frequencies. 

The dependence of cUp and cud on /3 for the viscoelastic 
Kelvin- Voigt model is given in Fig. [51 Here we used the 
analytical result for ujp, but LUd is determined numerically. 
This figure illustrates an important feature of the model: 
for P < I, ojd is relatively insensitive to /3, i.e. whether 
the system is solid or liquid, while Wp decreases linearly 
with (3. Thus, as /3 ^ 0, ujp/ojd -^ 0. This is consistent 
with the fact that as the system becomes more solid-like, 
the initial deviations from nonlinearity are positive, so 
the transition to liquid-like behavior is delayed to higher 
frequencies. We expect similar behavior for more general 
models. By measuring these characteristic frequencies, 
future experimental studies will be able to characterize 
the material properties of foams and other complex fluids. 

It is helpful to summarize the three characteristic 
frequencies — Wp, w„i, and Ud — that were defined above. 
LOp is the crossover frequency that characterizes the 
change in the scaling behavior of the total phase shift as a 
function of frequency. For pure solids, ujp is not defined, 
for pure liquids, cUp = 0, and for viscoelastic materials 
LOp > Q. ujni is the frequency above which we observe 
deviations from linear behavior in the horizontal veloc- 
ity profile. For uj > w„;, pure liquids display decaying 
nonlinear velocity profiles with positive curvature, and 
decay more strongly with increasing frequency. In solids 
and viscoelastic fluids, when oj > w^j, the horizontal ve- 
locity profile initially possesses negative curvature with 
deviations 'above' the linear velocity profile. Thus, the 
curvature of the profile at low frequencies near ujni can be 
used to differentiate 'liquid-like' from 'solid-like' velocity 
profiles. In the experiments, the solid-like response of 
the system at low frequencies is weak. Thus, we focus 
on measuring w^, the frequency above which the system 
begins to display decaying 'liquid-like' velocity profiles 
(instead of LOni) and tOp in the simulations and experi- 
ments. 



III. SIMULATION METHODS 

We performed numerical simulations in 2D of the bub- 
ble model, which was generalized to include particles with 
nonzero mass m, undergoing boundary-driven, oscilla- 
tory planar shear flow. The systems were composed of a 
total of Nt — 1280 disks with the same mass. Half of the 
disks were small and the other half were large with diam- 
eter ratio r = 2 to avoid crystallization under shear and 
match the bubble distribution used in the bubble raft ex- 
periments. The original simulation cell was rectangular 
with ^-coordinates in the range [0, L^] and y-coordinates 




FIG. 5: Plot of Up (solid line) and ujd (stars) versus /? = 
GpL^/rf for the viscoelastic Kelvin- Voigt model. Note that 
uop < LUd over a wide range of /3. 



in the range [— Ly/8, 9Lj,/8]. Bubbles were initialized 
with random initial positions within this rectangular do- 
main at a given packing fraction (p and then the system 
was relaxed to the nearest local potential energy mini- 
mum using conjugate gradient energy minimization [3l| 
with periodic boundary conditions in both the x— and 
J/— directions. All bubbles outside y = [0, Ly] formed two 
rough, rigid boundaries. Bubbles with y-coordinates in 
the range [Ly,9Ly/8] ([— Ly/8,0]) formed the top (bot- 
tom) boundary. N k, 1000 disks flUed the interior of the 
cell between the two rigid boundaries. After the top and 
bottom boundaries were formed, we used periodic bound- 
ary conditions only in the x-direction. Packing fractions 
in the range = [0.85, 0.9] were investigated. 

In the bubble model, bubble i experiences two pair- 
wise forces from neighboring bubbles j that overlap i: 
1. repulsive linear spring forces that arise from bubble 
deformation 






(9) 



and 2. viscous damping forces proportional to the rela- 
tive velocity between bubbles that arise from dissipation 
between the fluid walls 



pv 



-b{vi -Vj), 



(10) 



where e sets the energy scale for elastic deformation, 
Cij — {(Ti + iyj)/2 is the average diameter, r^ is the 
center-to-center separation between bubbles i and j, 
fij — rij/rij is the unit vector that points from the cen- 
ter of bubble j to the center of bubble i, and b is the 
damping coefficient. Note that when bubbles i and j do 
not overlap, the pairwise forces FT. = F^- = 0. 

The ratio of the damping to inertial forces can be 
expressed via a dimensionless damping coefficient b* = 
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FIG. 6: Average shear stress "E^y plotted versus the applied 
shear rate 7 for a system with — 0.86 and fe* = 2 undergo- 
ing steady planar shear flow in 2D simulations of the bubble 
model. The crossover shear rate at which the system transi- 
tions from quasistatic to power-law behavior is 7^ ~ 10~*. 



ba/y^em, where a is the small bubble diameter. Un- 
derdamped (overdamped) systems are characterized by 
b* < b* {b* > b*), where b* = a/2 for linear spring inter- 
actions. We studied both under and overdamped systems 
in the range b* — [0.1, 3]. The units of energy, length, and 
time in the simulations are e, a, and a/ y^e/m, respec- 
tively. 

The time evolution of the position r^ and velocity Vi of 
an interior bubble i can be obtained by integrating the 
equation of motion 



m- 









Fl 



(11) 



We employed standard Gear predictor-corrector algo- 
rithms to numerically integrate Eq. 1111 for the positions 
and velocities of the interior bubbles [32] ■ This simple 
'discrete element', short-range model for 2D foams al- 
lows us to quickly and efficiently generate an ensemble 
of configurations with a given set of external boundary 
conditions. 

Oscillatory planar shear was imposed by rigidly moving 
all of the bubbles that comprise the top boundary in the 
x-direction as a function of time according to 



Xb{t) = AsnY[ijjt), 



(12) 



while bubbles that comprise the bottom boundary re- 
main stationary. A = 0.8cr and lo are the amplitude and 
angular frequency of the sinusoidal driving. At this am- 
plitude, we did not observe any Tl events [26| over the 
entire range of driving frequencies studied. 

We calculated several important physical quantities in 
the simulations, including the phase shift of bubble x- 
displacements relative to the motion of the boundary and 
the horizontal velocity profiles of the bubbles. The x- 
displacements and velocity profiles reached steady state 



after a few cycles; thus, we began measurements after 5 
cycles and continued for an additional 15 cycles to calcu- 
late averages. The height dependence of the phase shift 
and velocity profiles were measured by partitioning the 
simulation cell into equal-sized rectangular bins centered 
at y with height Ay « 2 large particle diameters, with y 
measured from the bottom stationary boundary. 

To calculate the local phase shift 5{y), we averaged 
the bubble x-displacement Ua:{y,t) relative to the initial 
position over all bubbles within the bin located at y. We 
then fit the average bubble x-displacement in each bin 
to Ux{y,t) r^ sin(a;i — 6{y)) to determine the local phase 
shift S{y). We measured S{y) at several times during a 
given cycle to verify that it was time- independent, and 
the bubble motion was periodic. To measure the average 
horizontal velocity profile Vx{y-, t) of the interior bubbles, 
we used a binning procedure identical to that employed 
to measure 5[y). 

An important characteristic time (or frequency) scale 
in the bubble model is the shear rate 7c at which the sys- 
tem transitions from quasistatic behavior at low shear 
rate (shear stress Yixy oc 7*^) to highly fluidized behavior 
at high shear rate (S^y oc 7", where a > 0) when the sys- 
tem is driven by steady planar shear [l5| . This frequency 
scale has also been measured in the bubble raft experi- 
ments, and thus Wc = 27r7c can be used to normalize the 
crossover frequencies obtained in experiments and sim- 
ulations. To simulate systems undergoing steady planar 
shear, we employed the same boundary-driven method as 
described above except Xb (i) = Ly'^t instead of Eq. [T^ 
The flow curve for a system with (f) — 0.86 and 6* = 2 
is shown in Fig. [6l where the virial expression including 
dissipative forces was used to calculate the shear stress 
[32l |. To determine 7c, we calculated the median of the 
data point at which the flow curve first deviates by more 
than 10% from the power-law behavior and the previous 
data point at higher shear rate. For the flow curve in 
Fig. [51 we estimate 7c ~ 2 x 10"^. 7c was determined for 
each value of b and 6. 



IV. EXPERIMENTAL METHODS 

The experimental setup to apply oscillatory planar 
shear to bubble rafts includes three components: a rect- 
angular basin, oscillating paddle, and imaging system. A 
schematic of the experimental setup is shown in Fig. [71 
The basin has dimensions 37 cm by 15 cm and was filled 
to a depth of 5 cm with a surfactant solution. A paddle 
was located in the middle of the basin, leaving a span 
of 9 cm between it and the opposite wall. As illustrated 
in the schematic in Fig. [3 the ends of the system are 
"open" in the following sense. The entire basin is filled 
with bubbles, and the paddle only spans the central por- 
tion of the system. Furthermore, only the central third of 
the bubbles in the region covered by the paddle are used 
in the data analysis, and thus edge effects are minimized. 
The paddle was driven by an M-drive 23 stepper motor 
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Stationary Wall with grooves 

FIG. 7: A schematic of the experimental setup that applied 
oscillatory planar shear to the bubble rafts. The U-joint high- 
lighted in the figure is used to convert the rotary motion of 
the driving motor into oscillatory, planar shear. 



located outside the basin. A rotor and universal (U-) 
joint were used to convert the axial drive of the motor 
into linear sinusoidal motion. The amplitude of oscilla- 
tion was varied by changing the contact point between 
the rotor and U-joint. 

The bubbles were confined between a movable paddle 
and a fixed wall. The bubbles were constrained to move 
with the paddle using metal tabs that extended one bub- 
ble diameter into the system. The tabs were spaced ap- 
proximately every 5 bubbles. The fixed wall consisted 
of a series of square indentations approximately the size 
of a bubble. This fixed the bubble velocities at the sta- 
tionary wall to zero. Because the first row of bubbles 
at each wall is interspersed with elements to hold it in 
place, slight distortions of the bubbles prevented accu- 
rate measurement of their positions. Therefore, in the 
experiments we defined the location of y = and y — Ly 
to be the boundary between the first and second rows of 
bubbles at the stationary wall and paddle, respectively, 
instead of the location of the physical boundaries. For 
a sinusoidally oscillating paddle that is initially undis- 
placed, this gives the boundary condition at y = Ly-. 
Xb{t) = Asm(ujt), where A = 0.8a and tu are the ampli- 
tude and frequency of the driving. At this amplitude, no 
Tl events were recorded over the entire range of uj. 

The bubble raft consisted of a bidisperse distribution of 
bubble sizes: a 4 to 1 ratio of 2.5±0.3 mm to 5.3±0.5 mm 
diameter bubbles, which corresponds to a diameter ratio 
r w 2.1. The solution composition was 80% deionized 
water, 15% glycerol, and 5% miracle bubble (Imperial 
Toy Corp), which is a commercially available surfactant. 
The bubbles were produced by passing compressed nitro- 
gen through the solution via a needle. The diameter of 
the needle and nitrogen pressure determine the final size 
of the bubbles. To create bidisperse bubble mixtures, we 
used two needles with different sizes at constant pressure. 
After the bubbles were produced, we stirred the solution 
with a glass rod to break-up large-scale crystalline do- 
mains so that only short-range order persisted. We tuned 
the pressure to prevent multiple layers of bubbles from 
building up in the z-direction. 

The definition of the gas (or liquid) area fraction for 
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FIG. 8: (Color online) Shear stress Yj^y versus the applied 
shear rate 7 for the same bubble raft system described in 
Fig. [T] undergoing steady planar shear. We estimate the 
crossover frequency 7c ~ 0.065s~^ above which the system 
transitions from the quasistatic to the power-law flow regime. 



a bubble raft is somewhat imprecise because the bubbles 
form 3D structures on the water surface, which makes it 
difficult to define the amount of fluid in the walls. Also, 
using definitions based on Tl events present challenges 
because true vertices do not exist in bubble rafts, and 
therefore, minimum edge lengths are not well-defined. 
This is in contrast to fully confined two-dimensional sys- 
tems, such as soap films, for which more precise defini- 
tions of liquid area fraction exist ^^ . Therefore, for bub- 
ble rafts, one typically reports an operational definition 
of the gas area fraction, 0, as the average cross-sectional 
area of bubbles that is visible in the images. This is typ- 
ically done by applying a fixed cutoff to the images to 
separate pixels inside and outside of the bubbles. One 
expects that this method will provide values for ^ that 
approximate the definition of the packing fraction used 
in the bubble model. Using this operational definition, 
we find — 0.86 ± 0.04 for the bubble raft experiments. 
The error is estimated based on a range of choices for 
reasonable cutoff values. As we will show, one advantage 
of our studies is that they can provide a method for cali- 
brating our definition of the gas area fraction for bubble 
rafts with the packing fraction for the bubble model. 

To visualize the bubble motions, a 210 x 240 pixel CCD 
camera was held above the basin. The fioor of the basin 
was constructed from glass, which allowed the entire bub- 
ble raft to be illuminated from below by an electrolumi- 
nescent film (manufactured by Luminous Film Inc.). A 
frame rate of 60 frames per second was sufficient to cap- 
ture the bubble motion. The raw experimental images 
were filtered and thresholded to demarcate the spatial 
location of each bubble. Further details of the image 
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FIG. 9: (Color online) The normalized average horizontal ve- 
locity Vx{y)/vx{Ly) (averaged over times when the boundary 
velocity is maximum) as a fiinction of distance y from the sta- 
tionary bottom wall for several driving frequencies ui in the 
bubble raft experiments. At frequencies u < Ud ^ 4.0 s~^, 
the profiles are roughly linear, while above this characteristic 
frequency, liquid-like nonlinear profiles are observed. 



processing are provided in Ref. [23, [2^, [SJI • 

To measure the velocity profiles, we divided the sys- 
tem in the y-direction into ~ 20 equal-sized rectangular 
bins between y = and y = Ly. The instantaneous ve- 
locities of bubbles were then calculated by subtracting 
the center of mass positions of the bubbles in consec- 
utive image frames for each bin. We employed a PIV 
procedure to measure the time evolution of the bubble 
velocities and thus calculate the phase shift relative to 
the driving. The horizontal component of the velocity of 
the interior bubbles can be parameterized by the rms ve- 
locity Vrms = a/(w^(j/)), frequency of oscillation w, and 
local phase shift S{y) with respect to the moving top 
boundary. The horizontal component of the velocity of 
the interior bubbles is therefore given by 



Vxiy,t) ^ Vmagiy) COs{ujt ~ S{y)), 



(13) 



where v^ 



mag 



= V2t 



which allows us to calculate the 



local phase shift 6{y) relative to the moving boundary. 
In Sec. FVl below, we will show results for the total phase 
shift AEE6{0)~6{Ly). 

In the oscillatory planar shear experimental setup, we 
are not able to directly measure the shear stress and thus 
the flow curve. However, the flow curve has been mea- 
sured previously for a bubble raft with similar parameters 
(shown in Fig. [8]) f23|. We employed the same method as 
in the simulations to determine the crossover shear rate 
7c « 0.065 s~^ at which the system transitions from the 
quasistatic to the power-law flow regime. 



FIG. 10: (Color online) Total phase difference A plotted as 
a function of the driving frequency uu for the bubble raft ex- 
periments (black squares). The two lowest frequency data 
points (red triangles) represent the maximum possible phase 
shift that can be measured given the resolution of our exper- 
iment. The solid red and dashed blue lines have slopes 1 and 
2, respectively. The intersection of these two lines provides a 
rough estimate for the crossover frequency, LOp « 1.9±0.9 s~^. 



V. EXPERIMENTAL AND SIMULATION 
RESULTS 

In Figs. [^ and [TUl we show the results from the exper- 
iments in which bubble rafts were subjected to low am- 
plitude oscillatory planar shear. Fig. [5] displays the hor- 
izontal velocity profiles Vx{y) (measured at times when 
the boundary velocity is maximum as defined in Sec. lIVj) . 
over a range of driving frequencies. For this system, we 
find that the velocity profiles transition from nearly lin- 
ear to 'liquid-like' nonlinear near ujd ~ 4.0±0.5 s""'^. Note 
that at the lowest driving frequency, it appears that the 
velocity profile displays 'solid-like' behavior with slight 
negative curvature above a linear profile, but this feature 
is comparable to the size of the fluctuations. Since this 
feature is difficult to detect and the nonlinear liquid-like 
profiles are more robust in the experiments, we will focus 
on measurements of lu^ instead of ujni . 

In Fig. [ini we show measurements of the total phase 
difference A between bubbles near the top and bottom 
boundaries as a function of the driving frequency. As 
discussed in Sec.|lTl for a viscoelastic material we expect a 
crossover in the scaling of the total phase difference as the 
driving frequency is increased. For the experiments, we 
were unable to fully characterize the crossover behavior 
due to limits in our ability to measure the phase shift at 
low frequencies. However, using the maximum possible 
values of the total phase shift at the lowest frequencies 
red triangles in Fig. I10[ we can set a lower-limit on the 
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FIG. 11: (Color online) The normalized average velocity 
Vx{y)/vx{Ly) (averaged over times when the boundary ve- 
locity is maximum) as a function of distance y from the sta- 
tionary wall plotted for several driving frequencies ui in the 
bubble model simulations. For u) > uid ^ 0.02, the velocity 
profiles display 'liquid-like' nonlinear decay. 



exponent for the low-frequency scaling regime. This gives 
us a transition from A{lu) ~ w", with n > 2, at low- 
frequencies to linear scaling at high frequencies. Thus, 
for the bubble raft system, we estimate cOp « 1.9±0.9 s~^, 
and thus ujp < tOd- 

We also performed simulations of the bubble model 
undergoing small amplitude oscillatory planar shear over 
a range of packing fractions (j) and damping coefficients b* 
to compare to the bubble raft experiments. The trend for 
each set of </> and b* was similar: low frequency profiles 
are linear with a nonzero total phase shift followed by a 
transition to nonlinear liquid- like profiles when lo > LOd- 
In addition, the crossover in the scaling behavior of A 
versus lo was observed. This characteristic behavior is 
highlighted in Figs. [TT] and [T^ which show the velocity 
profiles and total phase shift for the bubble model with 
(/) = 0.86 and b* = 2. We find cod «^ 0.02 and tOp « 0.015 
for this set of parameters. In the low-frequency limit, we 
find A{u!) ^ Lo^, with n « 2.8 ± 0.3 (blue dashed line in 
Fig. [T^ , which is similar to the scaling predicted for the 
Kelvin- Voigt model |35j . 

Figure [T3] summarizes the simulation data (black 
squares) for the transition to liquid-like nonlinear profiles 
{ujd) and the crossover frequency (iOp). We show both 
characteristic frequencies uJd and cUp at fixed = 0.86 
as a function of b* including the under- and overdamped 
regimes and at fixed 6* = 2 as a function of 0. To directly 
compare the results in experiments and simulations, we 
normalize ujd and Up by Uc = 'injc obtained from the 
steady planar shear flow curves at each b and (j). The ex- 
perimental results are also displayed in Fig. [131 they are 
represented as solid lines since </> and b* are not known 
precisely for the bubble raft experiments. The simula- 
tions indicate that at fixed 0, the b* dependence of Wd/wc 
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FIG. 12: (color online) Total phase difference A plotted as 
a function of the driving frequency uu for the bubble model 
simulations at (f) — 0.86 and b* — 2. From the intersection 
of the solid red and dashed blue lines with slopes 1 and 3, 
respectively, we estimate cjp ^ 0.015. 



and LOpJLOc is fairly weak as the dissipation switches from 
underdamped to overdamped. At fixed &*, both lOdj^c 
and u>p/u>c increase with packing fraction, if we exclude 
the first point at </> = 0.85. 



VI. CONCLUSION 

We find that the bubble model captures all of the qual- 
itative features of the response of the bubble rafts to 
low-amplitude oscillatory planar shear flow. The key el- 
ements of the response include: 1. A regime with linear 
velocity profiles, but a nonzero total phase shift at low 
frequencies and 2. At the highest frequencies, a regime 
with liquid- like nonlinear profiles and nonzero total phase 
shift. In addition, in the low frequency regime, there is a 
crossover in the scaling of total phase shift as a function of 
the driving frequency. The bubble model reproduces each 
of these distinctive features found in the bubble rafts. 

The quantitative agreement shown in Fig. 1131 between 
the simulations and the experiments is also encouraging. 
In the range of parameters that we expect to correspond 
most closely to the experiments (0 « 0.86 and b* > V^), 
we find agreement to within error bars between the exper- 
iment and simulations for the characteristic frequency uJd 
at which the system transitions from linear to liquid-like 
nonlinear velocity profiles, and LOp, which characterizes 
the change in scaling of the total phase shift. 

The bubble model is often used to characterize highly- 
fiuidized flows where Tl bubble rearrangement events 
occur since it can quantify the statistics of these rear- 
rangement events. It has not been employed as often to 
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FIG. 13: (Color online) Summary of the measurements of the 
normalized frequencies uid/^c (top panels) and uip/coc (bot- 
tom panels) for the bubble model simulations (solid squares) 
and bubble raft experiments (solid lines), cod and ujp are nor- 
malized by u>c, which is the frequency at which steady planar 
shear flows crossover from highly-fluidized to quasistatic be- 
havior. For the simulations, we show ujdl^c and u)p/uc as a 
function of 6* at fixed (f> — 0.86 (left panels) and as a function 
of (ji at fixed b* = 2 (right panels) . Because it is difficult to de- 
fine (j) and b* precisely in experiments, the results for the bub- 
ble rafts are presented as solid lines (in between dashed lines, 
which give error bars for the experimental measurements) to 
identify the range of (j) and b* that best fit the experiments. 



quantitatively study slow, dense flows where bubble re- 
arrangements are rare since the results can depend on 
the dissipation model. One of the unique aspects of this 
study is that it allows a direct comparison between the 
bubble model and bubble raft experiments in the regime 
where Tl events and other large-scale bubble rearrange- 
ments do not occur, which enables stringent tests of the 
dissipation mechanism in the bubble model. The ob- 
served qualitative and quantitative agreement between 
simulation and experiment provides evidence that the 
bubble model provides a faithful, yet simple description 
of dissipation between liquid walls in bubble rafts. 

Our focus on the low-amplitude oscillatory shear 
regime, in the absence of Tl events, provides impor- 
tant insights into the response of foams to applied stress. 
First, we find that the bubble rafts clearly exhibit dissipa- 



tive behavior despite the absence of Tl events. This can 
only be due to motion in the fluid fllms, and emphasizes 
that further quantitative studies of fllm-level dissipation 
mechanisms are crucial to understanding foam rheology 
[23 |. Second, studies of the response of 3D foams and 
emulsions to oscillatory shear suggest that there are im- 
portant differences between the two systems in the low- 
amplitude regime [36, 37, 38]. Thus, studying the low- 
amplitude regime will highlight key differences in the me- 
chanical response of a variety of soft glassy materials. 
Finally, in this article we presented novel results that 
demonstrate the existence of two additional time (or fre- 
quency) scales associated with the response of foam to 
applied shear: tOp and ujd (or ujni)- Thus, we expect that 
these characteristic time scales will also play a key role in 
determining the flow curve and velocity profiles for foams 
undergoing steady planar shear flow. 

A quantitative comparison between the bubble raft ex- 
periments and bubble model simulations also provides a 
method to calibrate the gas area-fraction of the bubble 
rafts, which is notoriously difficult to measure in foam 
experiments. Of the three frequently used quasi-two di- 
mensional experimental setups (bubble rafts, bubble rafts 
with a top glass plate, and bubbles confined between two 
solid surfaces), it is most difficult to define the area frac- 
tion in the bubble rafts. By combining the bubble model 
simulations and Surface Evolver computations |39l . |40|, 
it will be possible to define an effective area fraction 
that will allow direct comparison between bubble raft 
experiments and other quasi two-dimensional glassy and 
jammed systems confined to surfaces and thin films. In 
the future, we will perform additional experiments over 
a range of packing fractions, surfactants that yield bub- 
bles with varied elastic properties, and liquid viscosities 
to investigate whether the predictions of the simple vis- 
coelastic model in Sec. |n]for the G- and 77— dependence 
oi iVni and Up are valid for the bubble rafts. 
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